Efficient estimation of energy transfer efficiency in light-harvesting complexes 
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The fundamental physical mechanisms of energy transfer in photosynthetic complexes is not yet fully un- 
derstood. In particular, the degree of efficiency or sensitivity of these systems for energy transfer is not known 
given their non-perturbative and non-Markovian interactions with proteins backbone and surrounding photonic 
and phononic environments. One major problem in studying light-harvesting complexes has been the lack of an 
efficient method for simulation of their dynamics in biological environments. To this end, here we revisit the 
second-order time-convolution (TC2) master equation and examine its reliability beyond extreme Markovian 
and perturbative limits. In particular, we present a derivation of TC2 without making the usual weak system- 
bath coupling assumption. Using this equation, we explore the long time behaviour of exciton dynamics of 
Fenna-Matthews-Olson (FMO) portein complex. Moreover, we introduce a constructive error analysis to esti- 
mate the accuracy of TC2 equation in calculating energy transfer efficiency, exhibiting reliable performance for 
environments with weak and intermediate memory and strength. Furthermore, we numerically show that energy 
transfer efficiency is optimal and robust for the FMO protein complex of green sulphur bacteria with respect to 
variations in reorganization energy and bath correlation time-scales. 

PACS numbers: 



I. INTRODUCTION 



Over the past few decades, there have has been significant 
interests in monitoring and simulating excitonic energy trans- 
fer in molecular systems 1ITI-H2I1. Recently 2D electronic spec- 
troscopy demonstrated that the excitation energy transfer in 
photosynthetic complexes could involve long-lived quantum 
coherence Il3l - fl9"ll . These experimental observations have 
lead to vigorous theoretical efforts to study quantum coherent 
dynamics in light-harvesting complexes ll20H29ll and observa- 
tions of environment-assisted quantum transport 1I20I, |2ll l24ll . 
Moreover, various ways for partitioning the contribution of 
quantum coherence to the energy transfer efficiency (ETE) 
have been explored ll23ll25i l30l - l33ll . In spite of many years of 
theoretical studies in energy transfer, the role of quantum ef- 
fects in the biological performance of photosynthetic systems 
is not fully understood. In particular, it is not known whether 
it is necessary to include quantum effects to demonstrate the 
optimal efficiency of these systems, predict the outcomes of 
ultrafast spectroscopic experiments [34], and explain the evo- 
lutionary path of the photosynthesis complexes II35I1 . 

The major difficulty in studying such complex open quan- 
tum system dynamics arises from the lack of an efficient 
method for simulation under realistic conditions. In the rel- 
evant biological systems, the system-bath coupling strength 
and free Hamiltonian parameters have typically comparable 
strength. In such cases the popular methods developed for the 
extreme perturbative regimes of weak system or weak envi- 
ronment break down, such as Forster energy transfer or Red- 
field/Lindblad formalisms |36]. Moreover, the stochastic mas- 
ter equation approaches lead to an incomplete description of 
the dynamics, for instance the Haken-Strobel-Reineker model 
i37l l38ll . treating environment as a classical white noise, is 
unable to capture finite temperature limit and detailed bath 
spectral density [36]. Additionally, the bath has often enough 



memory that all the master equation methods based on Marko- 
vian assumption become inadequate. A hierarchy of coupled 
master equations has been recently presented by Ishizaki and 
Fleming 1 26, 39], based on earlier works of Kubo and Tan- 
imura [40-42], that provides a general benchmark for simula- 
tion of light-harvesting complexes in all non-perturbative and 
non-markovian regimes with arbitrary accuracy. However, 
these general methods inevitably involve significant computa- 
tional overhead with increasing the size of system, temporal 
coherence, and in the low temperature limit. Thus, a variety 
of alternative approaches for simulation of open quantum sys- 
tems and light-harvesting complexes have been proposed to 
capture certain aspects of non-Markovian effects, quantum co- 
herence and/or beyond second order perturbation corrections 

JUS 

In this work, we demonstrate that the second-order time- 
convolution master equation (TC2) can be applied for ap- 
proximate calculation of Energy Transfer Efficiency (ETE) 
of complex open quantum systems interacting with a Gaus- 
sion environment beyond extreme Markovian and perturba- 
tive limits. We first introduce a set of approximations to de- 
rive the TC2 without the usual weak system-bath coupling as- 
sumption. Our study is based on the earlier work of Cao on 
the generalized Bloch-Redfield master equation [43]. Simi- 
larly, a post-perturbative derivation of Redfield equations can 
be obtained by employing harmonic oscillator baths and the 
Markov approximation 02811 . It has also been recently shown 
that, for dichotomic random two-jump process, a Stochastic 
Liouville equations can capture the excitonic dynamics be- 
yond the limits of both the weak coupling and short time cor- 
relation conditions |l58ll . A main task is to quantify how well 
such methods can capture dynamical properties of excitonic 
system in the intermediate regimes. Here, we provide a com- 
bination of analytical and phenomenological approaches for 
estimating errors, due to ignoring dynamical contribution of 



certain higher-order bath correlation functions, for calculation 
of ETE via TC2 master equation. This allows us to quantify 
the regions that TC2 would breaks downs for such compu- 
tation in the presence of environments with strong memory 
and strength. Our numerical simulations demonstrate that the 
estimated values of system-bath coupling strength and bath 
memory time-scale for the FMO protein indeed lead to opti- 
mal and robust energy transport in the intermediate regimes . 
This analysis could allow for a practical way to quantify the 
performance of large light-harvesting complexes and to ex- 
plore the optimality and robustness of these systems by rely- 
ing on a single time-convolution master equation. 

In a companion manuscript [60], we use this technique to 
comprehensively explore the efficiency of FMO protein and 
other small-size light-harvesting geometries. We demonstrate 
their optimality and robustness with respect to all the rele- 
vant internal and environmental parameters including, multi- 
chromophores spatial compactness, number of chromophores , 
spatial connectivity, dipole moments orientations, disorders, 
excitonic band gap structures, reorganization energy, tempera- 
ture, bath spatial and temporal correlations, initial excitations, 
and trapping mechanisms. We address whether or not the 
FMO complex structure and typically non-perturbative and 
non-Markovian environmental interactions are necessary for 
its performance. Specifically, we explore the general design 
principles for achieving optimal and robust excitonic energy 
transport and whether there are fundamental reasons for the 
convergence of time scales with resect to internal parameters, 
environmental interactions, and trapping mechanisms. 

This article is organized as follows. Section II discusses 
general energy transport dynamics and its modeling in com- 
plex quantum systems. In section III, the definition of ETE 
is introduced and for FMO complex, the optimality and ro- 
bustness of ETE as a function of system-bath coupling and 
bath memory is presented. In section IV, we evaluate the ac- 
curacy of our energy transport simulation. Detailed mathe- 
matical derivation of the TC2 equation and the error analysis 
for the approximate estimation of ETE are presented in the 
appendices. 



II. EXCITON TRANSPORT IN COMPLEX QUANTUM 
SYSTEMS 



Excitons are quasiparticles, each formed from a pair of 
electron and hole, that provide a natural means to convert 
energy between photons and electrons. There are two well- 
studied theoretical regimes of exciton transport. One ex- 
treme limit is semiclassical Forster theory in which the elec- 
tronic Coulomb interaction among different chromophores 
is treated perturbatively compared to the typically strong 
electron-phonon coupling. This condition leads to incoherent 
classical walks of exciton hopping among the chromophores. 
In the other extreme limit, the electron-phonon interaction is 
treated perturbatively using Redfield or Lindblad dynamical 
equations H20fl . However, the biologically relevant but less- 
studied cases are within the intermediate regime when the 
strength of the Coulomb and electron-phonon interactions are 




FIG. 1: The disordered structure of the Fenna-Matthews-Olson 
(FMO) complex: It is a trimer consisting of three identical monomers 
each formed by seven Bacteriochlorophylls (BChl) embedded in a 
scaffold protein. The FMO complex acts as an energy transfer 
channel in green sulphur bacteria guiding excitons from the light- 
harvesting antenna complex, in the proximity of BChls 1 and 6, to 
the reaction center which is in the proximity of BChls 3 and 4 



comparable. The HEOM in principle can be applied to all 
of regimes Il26i l36l l40i - l42ll . However, this approach is inef- 
ficient for exploring the properties of large light-harvesting 
complexes over a wide range of parameters, and for other ap- 
plications such as optimal material design IrSlll . 

Here, in order to arrive at a more efficient but less accurate 
approach, we start by considering the general time-evolution 
of open quantum systems. The dynamics of a photosynthetic 
system is influenced by the surrounding scaffold protein and 
solvent. Such an environment can be modeled as a phonon 
bath consisting of a set of harmonic oscillators. The total 
system-bath Hamiltonian can be written as 



total 
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where 
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Hs-ph 
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Here \j) represents the single exciton state of site j. The 
variables u>^, pj t % and qj£ are the frequency, position and 
momentum operators of the oscillator, respectively. The di- 
agonal elements {ej}s include system internal site energies 
plus reorganization energy shifts Xj = fioj^cfj induced 
by coupling to the phonon bath where d 3 is the dimension- 
less displacement of the (j, £)th phonon mode from its equi- 
librium configuration. The off-diagonal coefficients {Jj,k}$ 
represent dipole-dipole interaction between chromophores in 
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different sites. We assume that each site is linearly inter- 
acting with a separate bath, with operators Sj = and 
Bj = — fiugdj^Qj^ accordingly being the system and 
bath operators. A more general Hamiltonian should include 
system-bath coupling terms SijBij w i m the system op- 
erators Sij = however the fluctuations of the inter- 
dipole coupling amplitudes are typically much smaller than 
the site ener gy fl uctuations therefore we ignore these terms in 
our analysis lH, |63t] . 

The dynamics of an open system is given by quantum Li- 
ouvillian equation 



dp(t) 
Bt 



= £total[Ptotal(t)] — —ih([Ht tah Ptotal\t)])ph (3) 



with ptotai denoting the system-bath quantum state, and 
(...) p h being an average over phonon bath degrees of free- 
dom. The Liouvillian superoperator Ctotai is the sum of su- 
peroperators £5, C p h and Cs- P h corresponding to H$, H v h 
and Hs-ph- The explicit form of Ctotai can be obtained if 
the system and bath start in a product state: ps-ph(0) = 
p(0) <g> jOph(O). Furthermore, we assume that the phonon 
bath is initially in thermal equilibrium state at temperature T, 
exp(-PH ph )/Tr(exp(-/3H ph )) where /? = 1/kT. The as- 
sumption of an initial product state can be justified as the pho- 
tosystem is in its electronic ground state prior to interaction 
with a light source; if a quantum system is in a pure state, that 
implying a product system-bath state fl66ll . 

In the interaction picture, the compact formal solution of 
Eq. (O is 



P(t) = (7+ ex P / £totai(s)ds ) ph p(0) 



(4) 



where O denotes the interaction picture representation for an 
operator O. Expansion of the above time-ordered exponential 
function yields a Dyson expansion for time evolution of the 
density operator involving high-order bath correlation func- 
tions. Hereon we drop the subscript total and ph for simplic- 
ity. 



p(t)=y) f dh... f dt n (£(t)£(t 1 )...£(t n ))p{0) 

n J JO 

(5) 
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where the n-time correlation superoperator has the following 
form 

(£(t!)..x(t„))p(o) = H) B £ J2 C- 1 )""* 

ji..j„ ix.An 

x {B jh+1 (U k+1 )...B jn (UjB n (t h )...B jk (U k )) 

x S n (t n )...S Jk (t tk )p(0)S jk+1 (t ik+1 )...S jn (t in ) (6) 

with the site indices {ji,..-,j n } and the second summation 
over all indices {ii, i n } e {1, .., n} such that the t^ ti k 
and t ik+1 , t in are ordered backward and forward in time, 
respectively. The bath correlation functions vanish in a finite 
time interval and the system operator in each term of the above 



expansion Eq. (O has a bounded norm. Under these condi- 
tions the Dyson expansion always converges [67]. 

For the considered system-bath interaction (Q3 the bath op- 
erators B(tj) satisfies Gaussian statistics. That is, n-time cor- 
relation functions vanishes if n is odd. For even n, according 
to Gaussian property, the terms up to the second order in the 
bath correlation function are sufficient to exactly describe the 
dynamics of the system: 

(B(t H )...B(U 2n )) = £ H(l+B(t lk )B(t H )) (7) 

pairs l.k 

where the index pairs denotes the division of the labels 1 to 
2n into n unordered pairs. The operator I + is the index or- 
dering operator preserving the order of operators on RHS of 
Eq.fQ similar to the LHS. Note that here we have applied a 
generalized Wick's theorem in the form of Wightman func- 
tions rather than the usual form with time-ordered correla- 
tion functions ll64l l65ll . This expansion is possible assum- 
ing {Bj{ti)Bj{t 2 )) = (Bj{t 2 )Bj(ti}}* in consistence with 
Kubo-Martin-Schwinger condition [70]. The most general 
method to solve the master equation ([3]l is to utilize a path 
integral formalism leading to HEOM B42I1 . Here, we would 
like to avoid such a general approach since the required com- 
putational resources grows rapidly with the increasing size of 
the system as well as decreasing bath cut-off frequency and 
ambient temperature. 

In order to obtain a numerically efficient approach for sim- 
ulation of complex excitonic systems we first need to under- 
stand how the computational inefficiency arises. It should be 
noted that the source of computational complexity here is dif- 
ferent from that faced in quantum chemistry and condensed 
matter physics in ab — initio calculations of the ground state 
energy of interacting many-body fermionic systems. In such 
cases, the Hilbert space grows exponentially with the number 
of particles and the degrees of freedom involved. However, 
efficient techniques such as density functional theory B68I1 and 
density matrix renormalization group [69] can provide ap- 
proximate solutions. Here, due to the low intensity incident 
light and the time-scale separation of recombination process 
(in 1 ns) with fast energy transfer process (in 1 ps), we can 
ignore transitions between multi-excitation levels and focus 
on energy transport dynamics in a single exciton manifold. 
Thus, the Hilbert space grows linearly with the number of 
chromophores. However, due to the open nature of these com- 
plexes interacting with an environment that has strong mem- 
ory and strength, the time-nonlocal features of the dynamics 
are extremely difficult to simulate and explore. In these cases, 
reducing the number of environmental degrees of freedom 
(e.g, having a smaller bath frequency cutoff) does not translate 
into lower computational complexity. On the contrary it will 
enhance the non-Markovian behavior of the system. Specif- 
ically, the computational cost raises when attempting to treat 
these systems non-perturbatively while mapping the memory 
effects (e.g. encoded into a time-nonlocal kernel) to a set of 
coupled time-local master equations such as HEOM. In this 
work, we explore and quantify the applicability of TC2 mas- 
ter equation, as an inherently efficient method to capture cer- 
tain time non-local feature for calculation of ETE in environ- 
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ments with weak and intermediate strength. We first present 
a derivation of TC2 without the usual weak system-bath cou- 
pling assumption and quantify the errors in computing ETE as 
functions of both bath memory and strength. 

Here, we outline the main assumptions and steps of our 
derivation, for more details see Appendix A. First, we as- 
sume that the bath fluctuations are stationary; that is, these 
processes are insensitive to the reference point in time. For 
such quantum processes one can express bath correlations 
(Bj(t)Bj(t')) only as a function oft — t'. This character has 
also been assumed in HEOM. The correlation function is cal- 
culated from the bath spectral density 



Initial excitation: BChl 1 



Initial excitation: BChl 6 
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S(iU — U!^) 



{B 3 {t)B 3 {t')) = - 

7T 



dwJj(uj 



exp(—iuj(t — t')) 
1 — exp(— fihiu) 



(8) 



As we discuss in Appendix A, the stationary property of 
bath fluctuations can be exploited to provide a rather straight- 
forward solution for the equation of motion in the frequency 
domain. We obtain this time-nonlocal master equation by 
truncating the generalized Wick's expansion (0 to the sum of 
leading terms such that a two-point correlation can be factored 
out, that is 



(B(t il )...B(t i2 J) « (l+B(t 1 )B(t 2 ))(l + B(t k3 )..B(t k2 



(9) 
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FIG. 2: Evolution of site populations for all BChls of the FMO com- 
plex at T = 298K, X = 35cm -1 and 7" 1 = 166/ s. The results of 
simulations using HEOM published in I2H1 are shown in top panel 
(Courtesy of A. Ishizaki). Graphs A and B correspond to different 
initial states BChl 1 and BChl 6. The results of the simulation using 
TC2 are presented in graphs C and D are for initial states BChl 1 
and BChl 6. This comparison illustrates that the simulation by TC2 
yields oscillatory behavior as those simulations by HEOM approach, 
while significantly reducing the computational resources. However 
TC2 slightly overestimates amplitudes of oscillations in graph C and 
significantly underestimates the decay rates in graph D [28]. Never- 
theless, the measure of energy transfer efficiency that we employ in 
this work, as an important yield function to quantify the performance 
of light-harvesting systems (see section IV), is not very sensitive to 
actual oscillations. 



This approximation can be understood phenomenologically 
by noting that two-point correlation functions Cj(t,ti) = 
(Bj(t)Bj(ti)) typically decay with the time interval width, 
e.g., for a Dmde-Lorentzian spectral density, J(to) — 
2Xyu}/(uj 2 + 7 2 ), and at high temperature, they decay expo- 
nentially as e -7 ^ - * 1 '. The above approximation (0 is valid 
in the limit of large 7, however the resulting master equa- 
tion can capture some non-markovian behavior because of its 
time convolution form. Note that the relation (0 is differ- 
ent from the Bourret approximation for dichotomic process 
in which 2ri-point correlations can be exactl y expre ssed as 
(I?(fi)...i?(i 2n )) = IlLi<£(*2fc-i)£(*2fc)) Mlm Here, 
the remaining (2n — 2) point correlation (X + B(tk 3 )..B(tk 2rl )) 
in Eq. (0 still contains all the permutation of two-point cor- 
relation functions given by the expansion (J7). 

We can derive the TC2 equation by utilizing the above ap- 
proximations in the general equation of motion for the system 
density operator (0. In appendix A, we obtain the TC2 master 
equation 

j t p(t) = L sP {t) - Y}S 3 , i jf C 3 (t - t') x 

e Cs(t - t ">(S 3 p(t'))dt' - h.c], (10) 

where h.c. stands for Hermitian conjugate. The derivation 
of TC2 equation with Born (weak coupling) approximation is 



well known in the literature ll70ll . however, in the appendix A, 
we explicitly derive it under different assumptions away from 
weak system-bath limit. A phase-space representation for the 
above equation is also introduced in [43] as a generalization of 
the Bloch-Redfield equation which is equivalent to the second 
tier of hierarchy equations of motion (HEOM) in the high tem- 
perature and Lorentizan spectral density [39). Such general- 
ized Bloch-Redfiled equation is heuristically concluded from 
the Gaussian bath assumption which is inaccurate since it does 
not account for the pairing operation in the generalized Wick's 
expansion 0. Here, however, we constructively derive TC2 
by starting from the general Liouvillian equation and intro- 
ducing the required approximation (0 without making any 
weak ambient interactions assumption. Here we use TC2 to 
calculate ETE of a light harvesting complex equipped with 
an analysis to estimate the errors introduced by truncating the 
expansion (0, we find that the master equation can provide 
reliable estimations in the intermediate non-perturbative and 
non-Markovian regimes. An important feature of the TC2 is 
that it can be solved efficiently in time using Laplace trans- 
form technique. The numerical complexity of calculating the 
corresponding Laplace transform and its inverse does not de- 
pend on the form of the bath spectral density function or the 
low temperature limits in contrast to HEOM. 

We first use TC2 to simulate the quantum dynamics of 
the FMO complex (see Fig. at room temperature against 
HEOM as a general benchmark B26I1 . The FMO Hamiltonian 
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Initial excitation: BChl 1 Initial excitation: BChl 6 



No loss or trapping 



With loss and trapping 
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FIG. 3: Oscillatory dynamics of the FMO complex in the exciton 
basis at T = 298.K", A = 35cm" 1 and 7 _1 = 166/s. The energy 
eigenstates are spatially delocalized over various BChls, thus these 
oscillations manifest the presence of quantum dynamical coherence 
in the FMO complex. We observe that the exciton oscillations last 
for a few hundred femtoseconds endorsing experimental report of 
relatively long-lived quantum coherence beating. 
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is extracted from Ref.[62]. In both simulations the bath spec- 
tral density is considered to be Drude-Lorentzian with a two 
time correlation function of the form A(2//3 — i7)e~ 7 '* -tl ' 
where 7 and A are the cut-off frequency and reorganization en- 
ergy, respectively. Instead of solving over two million coupled 
differential equations needed for 1 1 hierarchy levels, as simu- 
lated by Ishizaki and Fleming [26], here we only need to solve 
the time-convolution equation ( fTOb in the frequency domain. 
Fig. |2] shows the oscillatory behavior of the population of all 
sites for two different initial conditions and a comparison with 
the results of Ref. l26ll . This data suggest that our simulation 
can capture the essential features of the FMO dynamics with 
a significant reduction in computational resources. It should 
be noted that our approach leads to slightly longer oscillations 
than those in HEOM. This is related to the known side effect 
of using TC2 yielding a double peak absorption spectra and 
it appears to be a generic artifact for any method relying on 
filtering or truncation of HEOM [HH Ell H ■ As we show be- 
low, this issue does not lead to a major problem in calculating 
energy transfer efficiency in which relies on the time average 
of the populations is involved. Similar simulations for cryo- 
genic temperature (T = 77K) are presented in appendix B. 

Our simulation of the dynamics of sites populations in Fig. 
[2] illustrated some oscillations that can last for a few hundred 
femtoseconds in agreement with the recent experimental re- 
sults at room temperature |fl9ll . However, population oscil- 
lations in the site basis per se cannot confirm the survival of 
quantum coherence. To this end, we simulate the dynamics 
in the exciton basis corresponding to the eigenvectors of the 
FMO Hamiltonian. We present the energy levels populations 
dynamics over time in Fig. ^indicating hundred of femtosec- 
onds oscillations. Next, we study the long-time dynamics of 
the FMO complex at the picoseconds time-scales and investi- 
gate its the thermal equilibrium properties. 



FIG. 4: Long-time dynamics of the FMO complex for initial excita- 
tion at BChl 1. Left panel presents the simulation without any lossy 
mechanisms. In contrast, the right panel illustrates the results in pres- 
ence of dissipation and trapping. Top/bottom panel are associated 
with dynamics in site and exciton basis. In graphs A and C, it can be 
observed that the FMO complex reaches to equilibrium state within 
lOps. In graphs (b) and (d) we note that the exciton is either fully 
absorbed or lost in the same time-scale of lOps. Thus, for most part 
the dynamics, the FMO complex is far from its equilibrium state. 

III. LONG-TIME BEHAVIOR OF EXCITONIC SYSTEMS 

Generally, there are two different competing electron-hole 
pair recombination processes that determine the energy trans- 
fer efficiency of light-harvesting systems. The first process 
happens within the time-scale of 1 ns due to dissipation (ra- 
diative and non-radiative) to the environment at each site. This 
adverse environmental effects guarantees that the ETE has a 
value less than one. The second recombination process is due 
to successful trapping at one or more reaction centers that typ- 
ically occur at the order of picoseconds. To capture the overall 
effect of electron-hole pair recombination at each site, an extra 
term C e - h is added to the TC2 dA"9V 

-p{t) = C s p{t) + C e - h p{t) (11) 
at 

-!>;< ^ Si C & ~ t')eZstt-t>) {Sjp{t , ))dtl _ hM] 

where C e - h = - £\ rj oss {\j){j\, .}-r tra p{\trap}(trap\, .}. 
Here \trap) is state of the site (BChl) connected to the reac- 
tion center. In this paper we assume that the reaction center 
is connected to the BChl 3 only: \trap) = |3). The coeffi- 
cients ri oss and r trap are the recombination and RC trap rates 
respectively and {, } is the anti-commutator symbol. Hereon 
we assume homogenous protein environments that all have 
similar correlation functions, Cj(t — t') — C(t — t'). Based 
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FIG. 5: Long-time dynamics of the FMO complex for initial excita- 
tion at BChl 6. Left panel presents the simulation without any lossy 
mechanisms. The right panel illustrates the results in presence of 
dissipation and trapping. Top/bottom panels are associated with dy- 
namics in site and exciton basis. Plots A and C illustrate that the 
FMO equilibrium is achieved within 5ps, faster than the equilibra- 
tion with initial excitation at BChl 1 as depicted in Fig. [4] This time 
difference can be understood by noting the presence of an energy bar- 
rier for excitation transport in the latter case. In graphs B and D the 
exciton is mostly absorbed or lost in less than 7ps, again implying 
non-equilibrium nature of the transport. 



on the dynamical equation ( fTTT i we provide a formal definition 
for energy transfer efficiency in the next section. 

Here, we study the equilibrium state of FMO dynamics us- 
ing Eq. (TTTT ). The long-time dynamics of an excitation ini- 
tially started at BChl 1 (6) is illustrated in Fig. g]© in both 
site and exciton basis. The left panel shows the dynamics in 
absence of any electron-hole recombination processes. The 
right panel includes the dynamics in presence of both dissipa- 
tion and trapping. In the absence of any lossy mechanisms the 
system relaxes to a equilibrium state p\ (oo) within lOps if the 
exciton starts from BChll (see Fig. |4](a) and (c)) and within 
hps if BChl 6 is the initial state (see Fig. [3] (a) and (c)). The 
longer equilibration time-scale of an initial excitation in BChl 
1 is a manifestation of an energy barrier on the exciton path to 
the trap site BChl 3. Such a barrier is missing on the exciton 
transfer from BChl 6 to BChl 3 path JH]. 

It is natural to assume that the FMO complex and vibra- 
tional mode of the scaffold protein are embedded in a thermal 
bath. Thus the combined pigment-protein complex should 
equilibrate to a thermal Gibbs state. Using TC2, we simu- 
late the infinite time behavior of FMO and find an equilib- 
rium state very close to its Gibbs state, denoted by po — 
exp(-PH s )/Tr(exp(-pH s )), for both initial state BChl 1 
and BChl 6: Trflp(oo) - pa\)/2 = 0.04. It should be noted 
that the true equilibrium state of FMO complex would be 



p G . = Tr ph (exp(-f3H to tai)/Tr(exp(-f3Htotai))- There- 
fore only for a very weak system-bath coupling the FMO 
steady state becomes pq. The TC2 captures this feature by 
noting that the distance of the equilibrium state from pc in- 
creases from 0.03 for A = 1cm -1 to 0.2 for A = 200cto~ 1 . 

If we include the electron-hole recombination processes 
due to loss and trapping, then as we expect the system relaxes 
to the (zero excitation) ground state. We observe that the re- 
laxation to the zero manifold occurs within lQps that is of the 
same order of time it takes for the system to equilibrate if the 
loss terms are ignored, see graphs (b) and (d) in Figs. |4]and|5] 
Thus, for the most part the processes of exciton energy trans- 
fer and trapping occur when the FMO complex is far from its 
equilibrium state. 



IV. ENERGY TRANSFER EFFICIENCY OF 
LIGHT-HARVESTING SYSTEMS 

A biologically relevant function for exploring the perfor- 
mance of light-harvesting complexes is the ETE as defined in 
Ref. Ir20l I22T f79tl . that is the total exciton population being 
successfully trapped. 



77 = 2r t 



rap 



(trap\p(t)\trap)dt 



(12) 



which is simply 2r trap {trap\p(s = 0)\trap) where p(s) is 
the Laplace transform of p(t). We provide a formal deriva- 
tion of the ETE in appendix C. At each moment, the over- 
lap of excitonic wave function with the site connected to the 
trap, (trap\p(t)\trap), quantifies the exciton availability for 
absorption by the reaction center. The ETE is indeed a sum- 
mation over probability of exciton presence weighted by the 
trapping rate, therefore a measure of successful transfer of the 
exciton captured by the reaction center. 

The method developed here allows us to efficiently sim- 
ulating the behavior of ETE as a function of various inde- 
pendent system and environmental degrees of freedom over 
a wide range of parameters. Two main parameters charac- 
terizing the effect of a Gaussian bath on an open system 
are the strength of the system-bath coupling and the bath 
internal memory time-scale. The former can be quantified 
by the reorganization energy shift Xj which is a quadratic 
function of the linear system-bath coupling strength given in 
Eqf2j The latter is determined by the width of the phonon 
modes spectral density (or cut-off frequency) denoted by 7,*. 
Here we assume that all Bchls are interacting with indepen- 
dent phonon baths with Drude-Lorentzian spectral density 
Jj(oj) = 2Xjjj0j/ (uj 2 + 7? ). Although this form of spectral 
density has been successfully employed for analyzing experi- 
mental results I72h74|1 . theoretical models suggest a summa- 
tion of Lorentzian terms with different A and 7 |T80L l8lll . We 
further assume that all baths have equivalent reorganization 
energy A and cut-off frequency 7. This assumption has been 
used in several empirical analyses ll62l |72T - T75ll . Here we ex- 
plore the variation of the ETE versus reorganization energy 
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Error in Calculating Energy Transfer Efficiency 



Bath Cutoff Freq. (cm 



FIG. 6: Energy transfer efficiency (ETE) of the Fenna-Matfhews- 
Olson Complex versus reorganization energy A (as a measure of de- 
coherence strength) and bath frequency cutoff 7 (as a measure of 
bath non-Markovianty). The bath frequency cutoff is plotted start- 
ing from 7 = 5cm" 1 due large errors of simulation in highly non- 
Markov regimes, see Fig. [7] The experimentally estimated values of 
T = 298K, A = 35cm" 1 , 7 = 50cm,- 1 or 150cm" 1 , r t -1 p = lps 
and r^ec = l ns reside at an optimal and robust neighborhood of 
ETE. For small reorganization energy (weak system-bath coupling 
strength) the non-Markovianity nature of the bath can slightly in- 
crease ETE. However, at larger reorganization energies it will sig- 
nificantly reduce ETE. Environment-assisted energy transport l24ll is 
clearly observed here. 



and bath frequency cutoff. Fig. [2] demonstrates the optimal - 
ity of ETE for the FMO protein complex at room tempera- 
ture T — 298K and at the experimentally estimated values 

of A = 35cm -1 , 7 = 50cm" 1 or 150cm -1 



assuming trap- 



ping and recombination rates as r^~ r = lps and r~ ec = Ins 
respectively consistent with Refs. iflll Ir32i f76l l79ll . It can be 
observed that the memory of the bath can slightly increase 
ETE in the regimes of weak system-bath coupling. Moreover, 
as expected the ETE drops significantly when the system in- 
teracts with a strong and slow bath. Here, the phenomenon 
of environment-assisted energy transport that was first sug- 
gested in the context of the Lindblad formalism (weak cou- 
pling and Markovian assumptions) [20] and Haken-Strobel- 
Reineker lf23l l37i l38tl . can be observed for all regimes us- 
ing our approach. An independent study on the optimality of 
ETE versus reorganization energy has been recently reported 
in Ref. [82]. The role of quantum coherence within the B800 
and B850 rings of LHII in purple bacteria for optimizing en- 
ergy transfer rates was also studied earlier by Jang, et. al. 
using generalized multichromophoric Forster theory Ref. ll46ll . 
We should mention that it takes about 1 .8 sec to calculate ETE 
using the method presented in this article on a desktop with a 
2.4 GHz processor and a 4 GB RAM memory. Overall, for 
a quantum system with the Hilbert space size d the computa- 
tional cost of simulations using TC2 grows as ad 12 for some 
constant a. 

Here we would like to discuss the connection between 
the reorganization energy and effective system-bath cou- 
pling strength. The reorganization energy, defined as 
J °° duiJ(u>) I (ttoj), is usually considered as the single param- 
eter characterizing the strength of the system-bath coupling. 
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FIG. 7: The error in calculating the energy transfer efficiency (ETE) 
of the FMO complex using TC2: this figure is a complement to Fig. 
[6]where ETE is plotted for different values of reorganization energy 
A and cut-off frequency 7. The error estimated by the function ( 117) 
increases with a smaller 7 and a larger A. Overall, the error anal- 
ysis demonstrates the reliability of the TC2 for estimating energy 
transport in excitonic systems beyond Markovian and perturbative 
regimes, specially for the estimated FMO values of A = 35cm -1 
and 7 = 50cm -1 or 150cm -1 marked by green lines. 



However, in a recent paper 115 911 . an insightful observation is 
made that this connection is not generally true and one should 
be careful when interpreting A as the magnitude of the sys- 
tem interaction with its environment. The reason is a sim- 
ple fact that the environmental oscillation modes far from 
the resonance frequencies of the system do not contribute 
to the decoherence dynamics of the system. Ref. 115 911 sug- 
gests to define an effective reorganization energy by A e / / — 
jvEmax ^j^j^ 1 (tjxj) where [E m i n , E max ] is the range of the 
relevant system frequencies. This range is [Ocm -1 , 550cm -1 ] 
for the FMO complex. We have examined this new quantity 
X e f / and considered the ETE as a function of A e / / and 7. As 
the simulation results presented in appendix E show, the dif- 
ference between ??(A e //,7) and r)(A,"f) becomes noticeable 
in the presence of a weakly coupled Markovian bath. 

In a companion paper (60J, we comprehensively explore the 
properties of the FMO complex as a function of all the rele- 
vant environmental and free Hamiltonian parameters. Next, 
we estimate the errors introduced into the calculation of ETE 
using the time-convolution master equation, and discuss the 
limitations and applicability of our approach. 



V. ERROR ESTIMATION OF ENERGY TRANSPORT 
SIMULATION 

The quantum dynamics as developed above was primarily 
motivated to lead to a time non-local master equation (TC2), 
Eq. < fm >. that may be readily solved in the frequency domain. 
We introduced a truncation of the correlation function expan- 
sion (0 by keeping the slowly decaying leading terms, and 
disregarding many fast decaying higher order bath correla- 
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tion functions, in order to arrive at a computationally efficient 
simulation of quantum dynamics via Laplace transformation. 
However, based on the above derivation, the regimes of the 
applicability of this method is not clear, since the errors intro- 
duced by such a truncation are not understood quantitatively. 
Although, it is qualitatively evident that ignoring such higher 
order bath correlation functions will introduce significant er- 
rors for very slow bath and very high reorganization energy. 
An important issue is the accuracy of this model in the in- 
termediate regimes. To address this issue one would ideally 
attempt to find an exact account of errors in various regimes 
of interest. However, this task is essentially entails calculat- 
ing the general evolution of the density operator of the sys- 
tem; that is the exact account of the errors is computation- 
ally as hard as simulating the the exact dynamics of system 
in all regimes. Nevertheless, using a combination of phe- 
nomenological and analytical approaches, we provide error 
estimation for weak and intermediate system-bath couplings 



and bath memory time-scales, thus quantifying reliability and 
applicability of our approach in these regimes. 

Here, we present an estimate of the ETE calculation inac- 
curacy associated to our approximation in Eq. ©. An upper 
bound for the error is 

POO 

Ai] = 2r trap \ (trap\p(t) - p T c2(t)\trap)dt\ (13) 
Jo 

where pit) and prc2 (t) is the exact density matrix of the sys- 
tem and ptci it) is the solution to the TC2. 

In order to estimate the above error we need to calculate 
pit) — prc2(t)- We use the Dyson expansion solutions for 
pit) and pTC2 it) in the interaction pictures in the absence of 
the term C e -h- The effect of this irreversible term is consid- 
ered later by introducing a decaying term exp(— r tra pt) (see 
appendix D). 



P~(t) ~ PTCl{t) = 

n J ° J ° h...j n n..i n 

(B jk+ AU k+ J..B hn it i2 jB h itiJ..^ (14) 

I 



where is the relative difference between coefficients 

of different terms in this expansion 

Aii. ..i„ _ i _ ( I + B h ( t h) B j2 {ti 2 )) (Z+B j3 {t l3 )..B ]2n jt l2n )) 

(Bj 1 (t il )...Bj an (t ian )) 

(15) 

The above expression for the difference of p(t) and p~TC2 it) 
is still exact. The coefficient A^'"? 2 ™ quantifies the relative 
error introduced into a 2n bath correlation function by keeping 
only the slow decaying leading terms within our approxima- 
tion. In order to compute the error in (Q~3}, we first need to 
estimate the error contributions from A; 1 '"?". 

Here we assume a Drude-Lorentzian bath, however the 
analysis can be simply repeated for other types of spectral 
density functions. Note that each term ignored in Q de- 
cays faster than the leading term. In appendix D, we exploit 
this feature to arrive at a computable form of A^"'| 2n inter- 
polating between time zero to infinity as 2 n-i (e~ 7 * + ... + 
e -(2n-i)-yty a j so use ^ following inequalityfTolto bound 
the time-ordered integral of bath correlation functions as 

| f dtx... f dt n {T+B(ti)...B(tan))ph\ < 
Jo Jo 

§^|2 f dtx I * dh(B{tx)B{ti))ph\ n (16) 

2 n n! Jo Jo 

where if^f is the number of contractions of 2n bath opera- 



tors. Appendix D shows how to take a further step to estimate 
the norm of system operators and arrive at an expression for 
an error estimation, defined by Eq. ( [T3V The estimate for the 
FMO complex is 

A V = mm(l, 2r trap £ / cfr ^ 1 e^"- 4 x 

1.75^.(1 + 4/( 7 /?) 2 )^ |* + -(e~ 7t - 1)1") (17) 

This estimation is for the high temperature limit 7 < See 
the appendix D for the error function in the limit of low tem- 
perature 7 > P^ 1 . The behavior of the above error function 
versus the reorganization energy A and cut-off frequency 7 is 
illustrated in Fig. [7] for a limited region of the ETE shown in 
Fig. [6] This figure shows that the TC2 equation can produce 
reliable results for intermediate values of the system-bath cou- 
pling and non-Markovian strength. The region in red denotes 
the parameter limits in which the TC2 can not produce a reli- 
able estimate of the ETE based on the error analysis. However 
as discussed in appendix D, the function ( fTTI i overestimates 
the error and a tighter bound may broaden the applicability 
regime of the TC2. The sharp transition from the blue region 
(almost zero) error to the red region (almost one) is due to 
the convergence/divergence properties of the series in func- 
tion ( fT7T >. Note that this error estimation is independent of 
the FMO complex properties. A rule of thumb for the appli- 
cability domain of the time nonlocal master equation can be 



9 



extracted from Fig. [7] for a given cut-off frequency 7, the 
system-bath coupling A should satisfy A < C7 where c = 6, 
for T — 298K and r^. ap = lps. For a detailed discussion on 
our error analysis see Appendix D. 

VI. CONCLUSION 

One of the main challenges in understanding the non- 
equilibrium behavior of photosynthetic complexes and de- 
signing artificial harvesting complexes is to efficiently sim- 
ulate their dynamics when the system-bath coupling is com- 
parable to the system energy scales. In this work, we have 
examined the TC2 master equation for simulation of exciton- 
ics dynamics as an efficient tool to compute ETE in complex 
quantum systems interacting with bosonic environments in 
low excitation limits. In particular, we have provided a deriva- 
tion of the TC2 master equation without making second-order 
perturbative assumption. We have applied this equation to 
calculate the energy transfer efficiency in Fenna-Matthews- 
Olson pigment-protein complex demonstrating optimality and 
robustness of energy transfer with respect to variations in bath 



temporal correlation and reorganization energy. In derivation 
of the dynamical equation some approximations have been 
made to allow a useful truncation of the n-time correlation 
expansion for bath time correlations. To account the inaccu- 
racies introduced by these approximations we have provided 
an error analysis estimating the parameter domain for applica- 
bility of TC2. Specifically, our study here quantify the errors 
in computing energy transfer efficiency of these systems via 
TC2. It will be of significant interest to similarly quantify the 
reliability of alternative schemes for calculations of popula- 
tion transfer such as Modified Redfield theory [83] or filtering 
of hierarchy of equations of motions HH. 
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Appendix A: Derivation of time nonlocal master equation 

The dynamics of a quantum system linearly coupled to a 
bosonic bath is generated by the quantum Liouvillian equation 



dp® 
dt 



-ih([H tot aup(t)])ph 



(Al) 



with Hamiltonian H given in Eq. (Q]). 

The interaction picture solution to Eq. (lAlt is 



/5(i) = (7^cxp / C to tai(s)ds ) ph p(0), 



(A2) 



where the interaction picture of any operator O is denoted by 
O. We consider the Dyson expansion of the above equation 
|p70h (Hereon we drop the subscripts "total" and "ph"): 



J 



P(t) 



dh / dt 2 {C{h)C{t 2 )) 



dt 1 



dU 



dU 



dU(C(h)C(t 2 )C(t 3 )C(t 4 )) 



p(0) (A3) 
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The time derivative of this expansion is: 



= 



dtx 



dt-, 



dh(C{t)C{t x )C{t 2 )C{tz)) 



P(0) 



(A4) 



r 



Note that for ti > t^.. > t n the n-time correlation superoper- 
ator has the following form 



(£(ti)...c(tn))p(o) = H)" E E 



i — k 



x (B jk+1 (Vi)«4 (*<»)% (*n)-% (*ij> 
x ^ )p(0)5, fc+1 ) ( A5 > 

with the second summation is over all indices {ii, ...,«„} £ 
{l,..,n} such that the U t > ... > t ik and U k+1 < ... < t in , 
and also with the site indices {ji, j„} 6 {1, .., N}. Here 
we assumed that the phonon bath is large enough to satisfy the 
Gaussian property; that is the 2m + 1-time correlation func- 
tions vanish and the 2m-time correlation functions can be de- 
termined by 2-time correlation functions: 

(B(t il )...B(t i2 J) = II (Z+B(ti)B(t k )) (A6) 

all l,kEii..i2m 
pairs 

where Z+ is the index ordering operator. This is a generalized 
Wick's theorem in the form of Wightman functions ll64l l65ll . 
We assume that the bath has a stationary memory function: 
{B j (t 1 )B j (t 2 )) = C j (t 1 -t 2 ). 
We keep the leading term in the above expansion I 



<B(t il )...B(t w )> « (l + B(t 1 )B(t 2 )){X + B(t k3 )..B(t k2 J) 

(A7) 

As it will be seen in appendix D, in this approximation we ig- 
nore the terms exponentially smaller than the remaining ones 
which is valid for relatively large values of 7. Note that the 
relation JA71 > is exact for a cross term due to the assumption 
of local baths, i.e. for j ^ j' 



(Bj )Bj {Ui)Bj> {U 3 )Bj> (Uj) - 



(A8) 



Now we can factor out f Q dt x {C{t)C{tx)) from Eq. (lA4l 
and obtain 



dt 2 



dt 3 {C(t 2 )C(h)) + ... P TC2(0) 



The operator pxc'2 represents the state of the system estimated 
by TC2. The above time non-local equation can be simply 
solved using the Laplace transform method. Thus, we have: 



-p^^W^^i))^^)^^),.] - h.c. 

3 

~^C j (t-t l )S j {t)[S j {t 1 ),.]-h.c. 

(A10) 



Now the explicit time convolution form of Eq. (IA9t in the 
Schrodinger picture becomes 



d_ 

at 



PTC2(t) = C S PTC2(t) - y^}Sj, —, 



h 2 



C 3 {t - t')e- iH ^ t - t ">/ h S J pTC2(t')e lHs{t - t ' )/h dt' - h.c' 



dt 1 (C{t)C(t 1 ))prc2{t\ 



(A9) 



(All) 



The above equation describes the influence of the bosonic 
(phononic) bath on the dynamics of the system. However, a 
photosynthesis complex is typically experiencing three other 
external dynamical processes, (a) The interaction with incom- 
ing light that generates the initial exciton in the system, (b) 
Recombination of the electron-hole pair at each site (BChl) 
leading to dissipation of the exciton to environment, (c) The 
trapping mechanisms that capture the excitation energy for 
charge separation in the reaction center which eventually con- 
verted to biochemical energy. In studying the excitonic dy- 
namics of a FMO complex the first process is negligible, since 
each FMO monomers has a very small absorption cross sec- 
tion and acts merely as an independent wire to transfer sun- 
light energy that is already absorbed by the antenna complex. 
The two other irreversible processes are modeled by adding 
two corresponding lossy terms to the RHS of Eq. (IA1 11 1 as: 
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§- t PTC2(t) = £ SP TC2(t) + £e- hP TC2(t) - £ [Sj , ^ J Cj(t - H) e -™s(t-f)/H g. p^^e^^^dt' - h.C 



(A12) 



where C e -h = - £\ r^ rec {\j)(j\, .} - r trap {\trap)(trap\, .} 
with r rec (rtrap) being the recombination (RC trap) rate and 
\trap) represents the state of the BChl connected to the RC. In 
this article we use Eq. (IA12b as the main dynamical equation 
describing the excitonic energy transfer process in the FMO 
complex. 

It should be noted that the well-known derivation of TC2 
is based on the assumption that system-bath coupling is so 
weak such that the third and higher order terms of the expan- 
sion ( IA3b and ( IA4b can be ignored. Basically, one can solve 
( IA3b to find p(0) as a function of p(t\) and substitute it in Eq. 
(1A4I >. all up to the second orders of C, and obtain TC2. In this 
paper we avoid the assumption of the weak system-bath cou- 
pling and we keep all powers of the reorganization energy in 
(1A41 >. Instead, here, we introduce the approximation (IA7t to 
arrive at Eq. dAl 11 1. Notice that in our approach the expanded 
solution of the density matrix has not been truncated over a 
finite power series of a "small" physical parameter (e.g. site 
energies, dipole-dipole interactions, reorganization energy, or 
system-bath Hamiltonain parameters in any fixed or rotating 
reference frame). 

Finally we would like to clarify that the approximation (|A7) 
is different from the Bourret approximation for which the 
Eq. (IA7t is applied recursively that is the (2n — 2) corre- 
lation function on RHS is also approximated 158, 17111 . Note 
that in contrast we consider a Gaussian bath and apply the 
approximation (lA7b separately to each term of the expan- 
sion (IA4b and thus treating the (2n — 2)-point correlation 
functions on RHS of Eq. JA7b exactly. That is in our ap- 
proximation we keep one out of the overall 2n terms in the 
Wick's expansion for the final time step. On the other hand, 
within the Bourret approximation, one attempts to capture the 
dynamics by keeping only one of the (2n/e)™ for 

large n) terms for two-point correlation functions. Neverthe- 
less, the Bourret approximation become exact when the noise 
can be regarded as a dichotomic random process leading to 

(B(ti)...B(t2n)) P h = Y[k=i(B(t 2 k-i)B{t 2 k)) P h- It remains 
an open problem to express the dynamics of the bath as a phys- 
ically motivated random process leading to the relation (IA7b . 
For a detailed error analysis of our approximation for comput- 
ing energy transfer efficiency see appendix D. 



Appendix B: Excitonic Dynamics at Cryogenic Temperature 

We simulated the excitonic dynamics of FMO at T = 77 K 
using TC2. Notice that here we consider 7 -1 = 50/sandA = 
35cm -1 to make a comparison with simulations in Ref. l26ll . 
At this temperature the bath two-time correlation function can 



Initial excitation: BCtil 1 Initial excitation: SChl 6 




a 200 400 600 800 a 200 -400 600 800 1000 




200 400 600 800 2O0 400 6BO BOO 1 QOO 



Time t (fs) Time f (fs) 

FIG. 8: Evolution of site populations for all BChls of the FMO com- 
plex at T = 77 K, A = 35cm -1 and 7 -1 = 50 fa. The results of 
simulations using HEOM published in 1 26] are shown in top panel 
(Courtesy of A. Ishizaki). Graphs A and B correspond to different 
initial states BChl 1 and BChl 6. The results of the simulation using 
TC2 are presented in graphs C and D are for initial states BChl 1 and 
BChl 6. TC2 predicts faster damping with initialization in BChl 1 
and almost no oscillation for BChl 6 as the initial state. 

be well approximated by C(t — t\) = A(i-(1 — ^1^2 ) — 
i 7 ) e -7(t-*0 + p{ ^. 2) 6(t), where v = §| Fig. © 
illustrates the simulation of BChls population for both initial 
state BChl 1 and BChl 6. Notice that in Fig. (O, 7 is almost 
three times larger than in the Fig. (|2), approaching Markovian 
regime. In comparison to HEOM results, TC2 underestimates 
the oscillations for initial excitation at BChl 1 and predicates 
negligible oscillations for initial excitation at BChl 6. 



Appendix C: Derivation of Energy Transfer Efficiency Function 

Generally, the sunlight energy captured by the antenna 
complexes has to be transferred to one or more reaction cen- 
ters for storage. However, there is always some finite chance 
of radiative or nonradiative electron-hole recombination at 
each BChl sites leading to dissipation of energy to environ- 
ment as florescence or quenching. The amount of the initial 
exciton that eventually arrives at RC determines the efficiency 
of the energy transfer process. 

Equation dA12b captures the dynamical evolution of the 
system in the single-exciton manifold. However, for a com- 
plete picture we consider the dynamical equation for the full 
quantum state Pf u u over an extended Hilbert space. This 
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(N + 2) -dimensional extended Hilbert space is constructed 
from N energy levels in single excitation manifold, the single 
state, |0), in zero excitation manifold, and finally a single state 
\RC) representing the reaction center. Modeling the presence 



of the reaction center by just an additional state is allowed 
since only exciton population and not coherence is transferred 
from the system to the reaction center. 



§- t Pfuu(t) = -i£sPfull(t) + £t"p f uu(t) - ^ J Cj(t - t l )e-^ t - t '^ h S ] p fu u{t l )e^ t -^' h dt' - h.c] 

(CI) 



where 



^trap 



{\trap){trap\,p} + 2\RC)(trap\p\trap){RC\ 

(C2) 



The total amount of the initial exciton population finally 
trapped in the reaction center is (RC\pf u u(oo)\RC). This 
can be equivalently evaluated by using Eqs. ( ICU and dC2b as 
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dt 

therefore 



(RC\p fuU (t)\RC) = 2r trap (trap\p fuU {t)\trap) (C3) 



7] = (RC\p fuU (oo)\RC) = 2r trap / dt(trap\p fuU (t)\trap) 



= 2r trap I dt{trap\p T C2(t)\trap) 
Jo 

(C4) 

Note that the RHS of the above equation is equivalent to 
2r tr ap(trap\p(s = 0)\trap), where p(s) is the Laplace trans- 
form of p(t). The convolutional form of the Eq. ( IA1 2b allows 
simple calculation of the p(s), which yields a closed form for 
the energy transfer efficiency function. 

The Laplace transformed form of the main equation (1A1 U 

is 

sp(s)-p(0) = C s p(s) + Ce-hP(s) 

J2[Sj,K j (s)-K^(s)} (C5) 



where 



K*(s) = -C(s + i£ s )(Sjp(s)). 



We can find the efficiency function r\ by solving this equation 
forp(s = 0). 



Appendix D: Error analysis for simulation of transport 
efficiency 

In the previous sections, we showed that the ETE of an ex- 
citonic system can be calculated efficiently by combining the 
TC2 and Laplace transform technique. However, this compu- 
tational simplicity inherently induces errors associated with 
ignoring certain higher order bath correlation functions. 

In this section we estimate the error in calculating the effi- 
ciency function r\ ( IC4b - An upper bound for the error is de- 
fined as: 



A77 = 2r trap | / (trap\p{t) - p TC 2(t)\trap)dt\ 
Jo 

POO 

<2r trap \{trap\p{t) ~ p TC2 {t)\trap)\dt (Dl) 
Jo 

The exact dynamical equation for the system in the pres- 
ence of recombination and trapping is obtained from Eq. (|3) 
by adding the term C e _hp(t) 



dp® 
dt 



= £ e - h p{t) - ih([H tota i, p{t)}) 



(D2) 



An interaction picture solution to Eq. (IAU is given in (1A31 >. 
but the term C e -h,p{t) will make it difficult to find a com- 
pact solution to Eq. (ID2I) . This complexity arises from the 
non-invertibility of the loss propagator operator exp(£ e _/ji). 
Instead, we consider p(t) and pTC2(t) to be solutions to Eqs. 
(lAlt and (IA1 Ik respectively, which do not include the effect 
of loss due to recombination and trapping in the reaction cen- 
ter. We incorporate these effects heuristically by considering 
that both p(t) and prc2{t) decay as exp(— r tra pt)- The new 
function to evaluate is 



Arj = 2r tr 



dte- rt ™-° 1 



^\{trap\p{t) - p T C2(t)\trap)\ 

(D3) 



The TC2 (IA1 lb can be solved for f>TC2(t) yielding 
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PTC2(t) 



i+ [ dh [ 1 dt 2 (C(t 1 )£(t 2 )) 
Jo Jo 

dt x r dh r dt 3 r dt^cit^tih))^)^)) 

o Jo Jo Jo 



P(0) 



(D4) 



and the difference between the exact and approximate solutions is: 



P{t) - PTC2{t) = 

E / *x- 



/■*2„-l 

JO 



■J2r l jk + 1--Jk 
l?.nlk + l Ik 



dhn E E (-i) n+fc <:;; 

ji...j„ i\..i n 

( B jk+i ( t i k+ i)-- B 32r, ( t i 2n )B jl (t h )..B jk (t ik )) x S n (U 1 )...Sj h (U k )p(0)S jk+1 (t ik+1 )...S jn (t in ) 



where 



(1+B jl (t n )B n (t l2 ))(l + B j3 (U 3 )..B hri (U 2n )) 

(- B jl(*il)—- B j2n(*<2„)) 



(D5) 



(D6) 



Here we continue the analysis for a Dmde-Lorentzian bath 
at high-temperature /3 , with cutoff frequency 7 and reor- 
ganization energy A, which is the case for FMO complex: 

C(t - ti) = (B(t)B(h)) = A(2//3 ± i 7 )e-T(*-*i) (± sign 



is determined by the order of t and t\). The analysis can be 
repeated for other types of spectral density functions. 

Each term ignored in Eq. (1A51 > decays faster than the lead- 
ing term 



\{l+B{t 1 )B(t 2 )){l+B{t 3 )...B{t 2n ))\ > \(l + B(t 1 )B{t k ))(l + B(t 2 )...B(t k _ 1 )B(t k+1 )...B(t 2n ))\ (D7) 
This behavior is evident in time ordered correlation terms 

(B(t)...B(t 3 )) = C(t - h)C{t 2 - t 3 ) + C(t - t 2 )C(h - t 3 ) + C(t - t 3 )C(h - t 2 ) 

= A 2 (2//3-i 7 ) 2 e- 7( *- t3) [e 7( * 1 - t2) +2e" 7( * 1 -* 2) ] (D8) 



where we approximate e 7 (* 1_t2 ) +2e~ 7 ^ 1 ~* 2 ^ with e 7 ^ 1 ^* 2 ^ 
A similar calculation can be done for higher order terms. This 
approximation introduces larger errors for smaller cut-off fre- 
quencies 7. Thus, stronger non-Markovian characteristic of 
the bath results in higher inaccuracy. For short times t we 
can make an estimation by assuming that all the terms in the 
expansion (lD8b have the same magnitude: 

\C(t - h)C(t 2 - t 3 ) + C(t - t 2 )C(h - t 3 ) 

+ C(t - t 3 )C(h - t 2 )\ ~ 3\C(t - h)C(t 2 - t 3 )\ (D9) 

or for the 2n'th terms 

\C(t - h)C(t 2 ,t 3 , ...) + C(t - t 2 )C(t u t 3 , ...) 

+ C(t - t 3 )C(h,t 2 , ...)| ~ (2n - l)|C(f - h)C{t 2 ,t 3 , ...)| 

(D10) 

This implies that the coefficient A 2 ™ is renormalized as 
X 2n /(2n — 1). This approximation becomes exact for long 



I 

times. Therefore the coefficient error A goes from 2 "~ 2 at 

& 271 — 1 

short-times to the value of zero at long-times. Here we con- 
struct functions A interpolating between short and long times. 
To this end, we make the following assumption: 

\C(t — t k +l)C(ti, t k ,tk+ 2 , t 2 n-\)\ave 
< e 7 *|C(i — tk)C(ti, tk—l, tk+l, t 2n -i)\ave, 

(Dll) 

, for time i, and on average for different values of ti, t 2n -i, 
that yields 

\C(t — tk)C(t\, tk-i, tk+i, t 2n ^i)\ ave 

< e^'^lCit - h)C(t 2 ,t 3 , ...)U- (D12) 

Note that this assumption is additional to the Gaussian prop- 
erties of the bath, and it holds if higher order bath corre- 
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lation terms in the expansion (IA6l l decay as an exponential 
function of the bath cut-off frequency 7. Using this relation, 
we find a computable form for the error function A that is 



J 



2 ^ = l{e~ lt + ... + e - ( 2 ™ -2 ) 7 *) interpolating between fg^f at 



time zero to zero at time infinity. We substitute this expression 
into ( lD5t and find an estimate for the efficiency error 



,-(2n-2) 7 t 



dt- 



2n - 1 



_ e -r trap t x 



f 

A77 = 2r trap 2_ \ I 

n>l J ° 

f dh f ... / dt n T+ (B(U k+1 )-B(U 2n )B(U 1 )...B(t lk ))\ 
\ zn )- Jo Jo JO 4 ■ 



\(trap\eM-iHst)T + S h (tiJ...S jh (^ (D13) 

I 



J1—32 



Here we used the identity: f Q dt\... f Q n 1 dt n C(ti)...C(t n ) — 

■^f J* dt\... L dt n T+C(ti)...£(t n ). The time ordering opera- 
tor, 7+, is not very suitable for our purpose here. This oper- 
ator is ordering the superoperators C(t) and not the bath or 



system operators B(t) or S(t). Instead the time ordering of 
C(t)s imposes a more complicated ordering of B(t) or S(t)s 
as described in the expansion ( |A5t . 

Next we consider an estimate for 



I (trap\ expi-iHs^T+Sj, (U^.-.S^ (t ik )p(0)S jk+1 (U k+1 )...S hn (t i2n ) exp(iH s t)\trap)\. 



31 ---Jr, 



(D14) 



Each operator S(tk) in the interaction picture equals 
cxp(—iHstk)S exp(iHstk)- Estimating this bound is not 
tractable in general, instead here we make an intuitive esti- 
mate. The coefficient Z is a product of the terms L(j,j') — 

\(j\eM-iHst k )\j')\ = VJ2 a \(M a )\ 2 \(j'\^a)\ 2 where 
\ip a )s are exciton basis. This expression represents how de- 
localized a site state can become over time. For a system with 
small terms L(j,j' ^ j) we can ignore the cross terms in 
( ID14t . This is the case for FMO with average L(j,j' ^ j) = 
0.05 B62I1 . From the average value of the terms L(j,j) = L 
we estimate Z ~ L 2n+2 . We consider L = 0.5 for FMO de- 
noting that the wavefunction on average becomes delocalized 
over two sites. In addition note that there are 2 2n different 
indices {ii..i n } and TV number of sites. 

By applying the Gaussian property (IA6t one can find rf83fl 



dti / ... / dt n (T+B(t 1 )...B(t 2n ))\ 







Jo 



< 



(2tQ! 
2 n n\ 



dt x / dt 2 {B{t 1 )B(t 2 ))\ r 



(2n)! 



where is the number of contractions of 2n operators 

B{h) to B(t 2n ). 

We now substitutes Z and ( ID15I ) into Eq. ( ID 13b to find the 



|A V '4/(7/3) 2 + l(i+-(e- 7 '-l))r i . (D15) 



final expression for the error estimate 



A77 = 2r trap ^2 



n>l 



dt 



E2n-2 - m -yt 
m=l e 

2n-l 



e ~r trap t x 



\n 1 

N4 n L 2n+2 — (1 + 4/( 7 /3) 2 )t It + ±( e -T* - l)|" (D16) 
n! 7 



Appendix E: Effective Reorganization Energy 



In a recent study H59H . Ritschel et al. introduce an effective 
reorganization energy as a measure for the system-bath cou- 
pling strength following the fact that the bath modes which 
are off-resonance with the system frequencies have negligi- 
ble effect on the dynamics of the system therefore the defini- 
tion of the reorganization energy A e / / = / °° duij(ui) / (irui) 

should be modified as X e ff = J^ mam dujJ{io)/{Tiuj), where 
E m ax and E m i n are the maximum and minimum system 
frequencies. For the Lorentzian spectral density, J(lo) = 
2Xjuj/(uj 2 +7 2 ), and [E min ,E max ] = [0cto~\ 550cm _1 ] 
we find A e // = (2A/7r)tan _1 (550/7). It can be seen from 
the notion of an effective reorganization energy that the bath 
modes and therefore a cut-off frequency 7 beyond 550cm -1 
are less relevant for the FMO excitonic dynamics. We have 
plotted the ETE and the estimated error for fixed effective re- 
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FIG. 9: The ETE as function of the effective reorganization energy. The sold (dashed) line shows ETE for a fixed value of A e / f (A) for three 
different values of 30, 60 and 100 cm -1 . The difference between ?7(A e //,7) and r](X, 7) becomes noticeable in the regime of large 7 and 
small A. The insets show the error estimation versus bath cutoff frequency for the fixed effective reorganization energy. 



organization energies in Fig. (0. A one percent difference be- 
tween T)(\ e ff, 7) and ?/(A, 7) plots is observed in the regime 
of large 7 and small A, as was expected. 

The energy transfer efficiency by definition (1C4) has a max- 
imum value of one. Therefore, any erroneous estimation of 
the efficiency should be less than one. However, in the above 
analysis the error is overestimated due to summing the magni- 
tude of all the terms. That might cause errors larger than one 
in some regimes which we reexpress them with a maximum 
error value of one. A similar approach is taken to estimate the 
noise threshold for fault-tolerant quantum computation in the 
presence of Gaussian noise 11851 18611 . where the threshold is 
found based on bath two-time correlation functions. 

It should be noted that the above analysis has been done 
for high temperature limit characterized by 7 < /3" 1 . In the 
landscape study presented in this paper, we also consider 7 



values greater than /3 _1 for which some correction terms need 
to be added to the bath correlation function. For T = 298° K 
and 7 < 500 cmT 1 , it is enough to consider the following 
expression as the bath correlation function [26, 87t 

(B(t)B(t% h = \d - ^'J/ - - z 7 )e-^-*') 



v /3 (2tt//3) 2 -7 2 



(2^//3) 2 -7 : 



;S(t-t') 



(El) 



Including this new correlation function, the error estimation 
can be simply obtained by replacing the coefficient (1 + 
4/( 7 /3) 2 ) with (1 + - (2 Jl% ir ?) inEq. (EH). 



